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Abstract 

A modification of the implicit algorithm for particle-in-cell simulations pro- 
posed by Petrov and Davis pQ is presented. The original lattice arrangement 
is not inherently divergence-free, possibly leading to unphysical results. This 
arrangement is replaced by a staggered mesh resulting in a reduction of the 
divergence of the magnetic field by several orders of magnitude. 
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1. Introduction 

In order to correctly reproduce physical processes in a particle-in-cell 
code, Maxwell's equations need to be solved consistently. However, the re- 
quirement of the magnetic field being divergence free is often violated by 
numerical algorithms leading to unphysical results [5J. Consequently, several 
divergence-cleaning schemes have been proposed, providing a way to remove 
magnetic source terms after the fact. Another possibility to correctly incor- 
porate Gauss' law for magnetism into a PiC-code is to use the staggered mesh 
first proposed by Yee j3] in 1966. The special arrangement of electric and 
magnetic fields inherently conserves a zero- valued divergence [I], provided 
that V--B = 0att = 0. A comparison by Balsara and Kim [5] identifies 
several problems of divergence-cleaning methods in MHD and notes their 
absence when using a staggered mesh. 
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Petrov and Davis pQ proposed an implicit particle-in-cell algorithm forgo- 
ing the staggered mesh approach. Continuing previous work by Kilian et al. 
[6] we intend to use this algorithm to study particle acceleration in astrophys- 
ical plasmas while keeping unphysical effects to a minimum. In this paper 
we therefore modify the scheme, incorporating the Yee lattice and effecting 
a reduction of V ■ B by several orders of magnitude. 



2. Definitions 

The quantities from pQ that are relevant to this paper are 

on+1/2 _ n a<lci rpn+1/2 /i\ 

and 



„ n+l/2 « 



1/2 = Uaq : +1/2 (f a + t n+ 1/2 (k x A^r 1/2 ) ) . (2) 

The tensor T (with indices suppressed for brevity) is defined as 



t 



1 + | Aft | 
with 



1 + Aft x Aft x Aft y + Aft z Aft x Aft z -Aft y ' 
Aft x Aft y - Aft z 1 + Att 2 y Aft y Aft z + Aft x 
Aft x Aft z + Aft y Aft y Aft z -Aft x 1 + Aft2 



f? n+l/2 . , 
A Qn+l/2 = QgBa At 
« n+l/2 O 



(3) 
(4) 



The quantities the charge, mass, and number density of 

(computational) particle a. 7^ +1//2 is the particle's relativistic gamma factor 
and Ba +1 ^ 2 its local magnetic field at time n + l/2. is the momentum of 
particle a at time n. 

The deposition of S and Sj on the grid and the interpolation of E and B 
to the particle position is achieved via a standard weighting function. Our 
algorithm makes use of the triangular shaped cloud (TSC) scheme. 



3. The modified lattice arrangement 

The original algorithm by Petrov and Davis [1] stores electric fields on 
grid nodes and magnetic fields in the cell center. The vector quantity Sj and 
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the tensor quantity S are deposited on grid nodes, as well. Since the electric 
field is updated according to 

(i + £ n+1/2 ) E n+1 =(l- S n+1 ^ E n +^ (V x H n+1 / 2 - Sj n+1 ^ , (5) 

and all required quantities are defined on grid nodes, this equation can be 
solved locally for E n+1 . 

Our approach keeps the original field layout by Yee [3] with the compo- 
nents of Sj stored like the corresponding components of the electric field. 5 
is stored on grid nodes and interpolated linearly for each component of the 
electric field to be calculated. When calculating the new value for e^ 1 ' 2 '^^ , 
S is taken to be (S^^ + S i+1 ' j ' k )/2, for Ej' j+1/2,k it is (S i ' j ' k + S i ' j+1 ' k )/2 
and for Ei' j ' k+1/2 it is {S^ + S i ' j ' k+1 )/2. 

Since S is not a diagonal tensor, all the components of Sj, V x B and 
E n need to be known at the same point as the component of E n+1 to be 
calculated, as well. These three quantities can be interpolated the same way. 

For £ x n+1 : 

A i+l/2,j,k _ £i+l/2,j,k /g-j 
£i+l/2,j,k _ (£i,j+l/2,k + £ij-l/2,k + ^i+l,j+l/2,k + £i+l,j-l/2,k^ ^ 
A i+l/2,j,k _ fj^i,j,k+l/2 _|_ £ij,k-l/2 _|_ ^i+l,j,fc+l/2 _|_ ^j+l,i,fc-l/2^ ^ ^ 

For £ y n+1 : 

= r^i+l/2,j,k _|_ ^i+l/2,j+l,fc + ^i-l/2J,fc _|_ ^i-l/2J+l,fc\ ^ ^ 

= A M+l/2,fc (1Q) 

£i,j+l/2,k _ /£ij,k+l/2 + ^i,j+l,fc+l/2 + ^lJ,fc-l/2 + ^i,j+l,fc-l/2^ ^ ^-g 



For £, 



n+l. 



^ij,fc+l/2 _ /^i+l/2J,fc + ^i+l/2J,fc+l + ^i-l/2j,fc + ^i-l/2,j,ft+l\ 
^i,i,fe+l/2 _ ^ij+l/2,fc + ^ij+l/2,fc+l + ^i,j-l/2,fc + ^i,j-l/2,fc+l^ ^ 
^i,i,fe+l/2 _ ^i,j,k+l/2 /-^ 

4. Simulation setup 

In order to compare the original lattice configuration with the Yee lattice, 
identical simulations are performed using both algorithms. 
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electron plasma frequency 
length of timestep 
Debye length 
cell edge length 
mass ratio 



CO- 



pc 



At 

A D 
Ax 

m p /m c 



2.0 • 10 8 rad/s 
4.1 ■ l(T 10 s 
7.5 cm 
21 cm 
42 



Table 1: Parameters of the thermal plasma. 



The size of the three-dimensional lattice is 64 x 64 x 64 with periodic 
boundary conditions. There are 25 electrons and 25 protons in each cell. Each 
simulation is running for 4000 timesteps. The plasma is chosen to be thermal, 
meaning that the velocity components for each particle are independent and 
follow a Gaussian distribution of width vth- The protons and electrons are 
in thermal equilibrium so that the width of the Gaussian distribution of the 
protons is equal to v t h divided by the square root of the mass ratio ^/m p /m e . 

The physical parameters are listed in table [TJ 



5. Results 

When directly comparing V ■ B at the end of the two simulation runs, a 
significant difference is manifest. The ratio of the total divergence of both 
simulations at timestep 4000 is 

EiifclV-Syeel _ 7 « 1n -L5 



Tliijk ' -^original | 



7.6 ■ 10~ 15 (15) 



Fourier transforming V ■ B at timestep 4000 in space and plotting the 
absolute values for a slice along the z-axis yields figure [TJ As can be seen, 
the original arrangement introduces high-amplitude patterns at short wave- 
lengths. The Yee-arrangement mostly produces low-amplitude noise, as de- 
sired. 

Furthermore, plots of the y-component of the magnetic field, prepared the 
same way, are shown in figure [2j The original lattice arrangement introduces 
unphysical behavior at short wavelengths which is not present when using a 
Yee-lattice. 

A dispersion plot of the y-component of the electric field along the x- 
axis for the simulation using the Yee-lattice is obtained as follows. For each 
timestep the E y values are summed up for slices perpendicular to the x-axis. 
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Figure 1: Spectral plot of the distribution of V • B in x- and y-direction for the original 
lattice arrangement (left) and the Yee arrangement (right). 




E y in x direction (transversal) 




A Fourier transformation in space and time results in a u-k x relation, de- 
picted in figure [3j In addition, the dispersion of the electromagnetic wave 
as predicted by theoretical calculations is shown (incorporating effects from 
finite At and Ax). As can be seen, the modified code correctly describes 
wave dispersion in a thermal plasma indicating that the alterations do not 
introduce unphysical effects. Likewise, further tests of more complicated sim- 
ulation setups do not show unexpected behavior resulting from our changes. 

Finally, the simulations are repeated with Aa; and At scaled down by a 
factor of 4, permitting a comparison with our explicit PiC Code ACRONYM 
[6]. The energy development for the implicit scheme using the Yee-lattice, 
the explicit scheme, and the unaltered implicit scheme are shown in figure 
|4j The original scheme shows deviations from the explicit results with the 
magnetic field energy growing steadily over the course of the simulation. In 
contrast, the altered scheme nicely reproduces the energy development of the 
ACRONYM simulation. 

As shown, it is possible to alter the algorithm proposed by Petrov and 
Davis pQ to obtain a divergence-free setup. In our code, the only steps that 
needed to be modified were the deposition and interpolation of grid quanti- 
ties and the calculation of the updated fields. The former change is easily 
implemented by adjusting the offsets used when evaluating the weighting 
function. The latter change is described in section |3j 

Due to the additional calculations introduced into the scheme, some per- 
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Figure 4: Energy development of the rescaled thermal simulation for the altered implicit 
scheme (top left), the explicit scheme used in ACRONYM (top right) and the original 
scheme (bottom). 

formance is lost. However, since most of the changes are confined to the field 
update step which only needs to be executed once per cell and timestep, the 
performance loss is small and on the order of a few percent. 
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